### Plot: Overview hybrid models

rm(list=ls(all=TRUE))


#packages
library(tidyverse)
library(sjmisc)
library(sjPlot)
library(sjlabelled)
library(ggplot2)
library(gghighlight)
library(RColorBrewer)
library(cowplot)
library(data.table)
library(dotwhisker)
library(ggpubr)
library(gridExtra)
library(lemon)
library(grid)
library(gginnards)

select <- dplyr::select
filter <- dplyr::filter


##### HYBRID MODELS ####

##### read data ####
setwd("SET_YOUR_PATH")

#hhinc_dec_pp_sqrt
bhps <- read_stata("./BHPS/bhps_coef_hhinc_dec_pp_sqrt_e.dta")
gles <- read_stata("./GLES/gles_coef_hhinc_dec_pp_sqrt_e.dta")
gss <- read_stata("./GSS/gss_coef_hhinc_dec_pp_sqrt_e.dta")
liss <- read_stata("./LISS/liss_coef_hhinc_dec_pp_sqrt_e.dta")
nlsy97 <- read_stata("./NLSY97/nlsy97_coef_hhinc_dec_pp_sqrt_e.dta")
polat <- read_stata("./POLAT/polat_coef_inc_dec_e.dta")
psid <- read_stata("./PSID/psid_coef_hhinc_dec_pp_sqrt_e.dta")
shp <- read_stata("./shp/shp_coef_hhinc_dec_pp_sqrt_e.dta")
soep <- read_stata("./SOEP/soep_coef_hhinc_dec_pp_sqrt_e.dta")

####p repare data ####

# delete empty rows
bhps <- bhps[-7,]
gles <- gles[-12,]
gss <- gss[-3,]
liss <- liss[-6,]
nlsy97 <- nlsy97[-3,]
polat <- polat[-10,]
psid <- psid[-2,]
shp <- shp[-4,]
soep <- soep[-6,]

#set new colnames 
colnames(bhps) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2", "id")
colnames(gles) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2", "id")
colnames(gss) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2", "id")
colnames(liss) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2", "id")
colnames(nlsy97) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2", "id")
colnames(polat) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2", "id")
colnames(psid) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2", "id")
colnames(shp) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2", "id")
colnames(soep) <- c("var", "estimate",  "std.error", "ci_lower", "ci_upper", "N", "r2", "id")


#add identifiers for datasets
bhps$df <- 1
gles$df <- 2
gss$df <- 3
liss$df <- 4
nlsy97$df <- 5
polat$df <- 6
psid$df <- 7
shp$df <- 8
soep$df <- 9




###recode DVs

#hhinc deciles (eq)
bhps$dv <- rec(bhps$id, rec = "1=3 [Political Interest]; 2=1 [Vote]; 4=2 [Intention to vote]; 
               3=4 [Internal Efficacy]; 5=5 [Duty to vote]; 6=6 [Vote: pers. benefits]")
bhps$term <- as_character(bhps$dv)

gles$dv <- rec(gles$id, rec = "1=1 [Vote]; 2=3 [Political Interest]; 3=2 [Intention to vote]; 
               4=4 [Has Party ID]; 5=5 [Party ID strength]; 6=10 [Internal Efficacy]; 7=11 [Imp.: Elections]; 
               8=12 [Imp.: Elec. Campaign]; 9=13 [Pol. media use]; 
               10=15 [Pol. Know.: Verbal]; 11=16 [Pol. Know: Visual]")
gles$term <- as_character(gles$dv)

gss$dv <- rec(gss$id, rec = "1=29 [GSS: Interest intl. affairs]; 2=30 [GSS: Interest mil. policy]")
gss$term <- as_character(gss$dv)

liss$dv <- rec(liss$id, rec = "1=3 [Political Interest]; 2=1 [Vote]; 3=2 [Intention to vote];  
               4=4 [Internal Efficacy]; 5=5 [Pol. part. index]")
liss$term <- as_character(liss$dv)

nlsy97$dv <- rec(nlsy97$id, rec = "2=2 [NLSY97: Political Interest]; 1=1 [NLSY97: Vote]")
nlsy97$term <- as_character(nlsy97$dv)

polat$dv <- rec(polat$id, rec = "3=2 [Political Interest]; 1=1 [Intention to Vote]; 2=7 [Pol. media use]; 4=3 [Internal Efficacy]; 
                5=6 [Duty to vote]; 6=4 [Has party ID]; 7=5 [Party ID strength]; 8=8 [Pol. part. index];  9=9 [Pol. part. onl.]")
polat$term <- as_character(polat$dv)

psid$dv <- rec(psid$id, rec = "1=1 [PSID: Vote]")
psid$term <- as_character(psid$dv)


shp$dv <- rec(shp$id, rec = "1=1 [Vote]; 2=2 [Political Interest]; 3=3 [No. polls part.]")
shp$term <- as_character(shp$dv)

soep$dv <- rec(soep$id, rec = "1=1 [Vote]; 2=3 [Political Interest]; 3=25 [Has party ID]; 
               4=2 [Intention to vote]; 5=26 [Party ID strength]")
soep$term <- as_character(soep$dv)


us <- bind_rows(gss, nlsy97, psid)
d <- bind_rows(bhps, gles, gss, liss, nlsy97, polat, psid, shp, soep)







#reorder rows for dwplot
bhps <- bhps[c(5,3,6,4,2,1),]
gles <- gles[c(11,9,10,8,7,6,5,4,3,2,1),]
liss <- liss[c(4,3,5,2,1),]
polat <- polat[c(9,7,4,3,6,8,2,1),]
shp <- shp[c(3,2,1),]
soep <- soep[c(5,2,4,3,1),]
us <- us[c(5,4,3,2,1),]



#multiply effects of LPM by 9 to make estimates comparable
bhps$estimate <- ifelse(bhps$dv == 1, bhps$estimate * 9, bhps$estimate)
bhps$std.error <- ifelse(bhps$dv == 1, bhps$std.error * 9, bhps$std.error)

gles$estimate <- ifelse(gles$dv == 1 | gles$dv==4, gles$estimate * 9, gles$estimate)
gles$std.error <- ifelse(gles$dv == 1 | gles$dv==4, gles$std.error * 9, gles$std.error)

liss$estimate <- ifelse(liss$dv == 1 | liss$dv == 3, liss$estimate * 9, liss$estimate)
liss$std.error <- ifelse(liss$dv == 1 | liss$dv == 3, liss$std.error * 9, liss$std.error)

polat$estimate <- ifelse(polat$dv == 4, polat$estimate * 9, polat$estimate)
polat$std.error <- ifelse(polat$dv == 4, polat$std.error * 9, polat$std.error)

shp$estimate <- ifelse(shp$dv == 1, shp$estimate * 9, shp$estimate)
shp$std.error <- ifelse(shp$dv == 1, shp$std.error * 9, shp$std.error)

soep$estimate <- ifelse(soep$dv == 1 | soep$dv == 25, soep$estimate * 9, soep$estimate)
soep$std.error <- ifelse(soep$dv == 1 | soep$dv == 25, soep$std.error * 9, soep$std.error)

us$estimate <- ifelse(us$dv == 1 | us$dv >= 22, us$estimate * 9, us$estimate)
us$std.error <- ifelse(us$dv == 1 | us$dv >= 22, us$std.error * 9, us$std.error)



#### create plots ####

#objective: hhinc_dec_pp_sqrt
p_bhps <- dwplot(bhps, vline = geom_vline(xintercept=c(0), color = "grey60"), dot_args = list(aes(color = "black")),
                 whisker_args = list(aes(color = "black"))) +
  scale_colour_grey(start = .1, end = .1)  + 
  geom_vline(xintercept=0.05, linetype = "dashed", color = "grey60") +
  geom_vline(xintercept=0.1, linetype = "twodash", color = "grey60") +
  scale_x_continuous(limits = c(-0.1,0.2), breaks = c(-0.1,-0.05,0,0.05,0.1,0.15,0.2), labels = c("-0.1","-0.05","0","0.05","0.1","0.15", "0.2")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("UK (BHPS & UKHLS)")  
p_bhps <- move_layers(p_bhps, "GeomPoint", position = "top")
p_bhps <- move_layers(p_bhps, "GeomErrorbar", position = "top")


p_gles <- dwplot(gles, vline = geom_vline(xintercept=c(0), color = "grey60")) +
  scale_colour_grey(start = .1, end = .1)  + 
  geom_vline(xintercept=0.05, linetype = "dashed", color = "grey60") +
  geom_vline(xintercept=0.1, linetype = "twodash", color = "grey60") +
  scale_x_continuous(limits = c(-0.1,0.2), breaks = c(-0.1,-0.05,0,0.05,0.1,0.15,0.2), labels = c("-0.1","-0.05","0","0.05","0.1","0.15", "0.2")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("Germany (GLES)")
p_gles <- move_layers(p_gles, "GeomPoint", position = "top")
p_gles <- move_layers(p_gles, "GeomErrorbar", position = "top")


p_liss <- dwplot(liss, vline = geom_vline(xintercept=c(0), color = "grey60")) +
  scale_colour_grey(start = .1, end = .1)  + 
  geom_vline(xintercept=0.05, linetype = "dashed", color = "grey60") +
  geom_vline(xintercept=0.1, linetype = "twodash", color = "grey60") +
  scale_x_continuous(limits = c(-0.1,0.2), breaks = c(-0.1,-0.05,0,0.05,0.1,0.15,0.2), labels = c("-0.1","-0.05","0","0.05","0.1","0.15", "0.2")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("Netherlands (LISS)")
p_liss <- move_layers(p_liss, "GeomPoint", position = "top")
p_liss <- move_layers(p_liss, "GeomErrorbar", position = "top")

p_polat <- dwplot(polat, vline = geom_vline(xintercept=c(0), color = "grey60")) +
  scale_colour_grey(start = .1, end = .1)  + 
  geom_vline(xintercept=0.05, linetype = "dashed", color = "grey60") +
  geom_vline(xintercept=0.1, linetype = "twodash", color = "grey60") +
  scale_x_continuous(limits = c(-0.1,0.2), breaks = c(-0.1,-0.05,0,0.05,0.1,0.15,0.2), labels = c("-0.1","-0.05","0","0.05","0.1","0.15", "0.2")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("Spain (POLAT)")
p_polat <- move_layers(p_polat, "GeomPoint", position = "top")
p_polat <- move_layers(p_polat, "GeomErrorbar", position = "top")
                 
p_shp <- dwplot(shp, vline = geom_vline(xintercept=c(0), color = "grey60")) +
  scale_colour_grey(start = .1, end = .1)  + 
  geom_vline(xintercept=0.05, linetype = "dashed", color = "grey60") +
  geom_vline(xintercept=0.1, linetype = "twodash", color = "grey60") +
  scale_x_continuous(limits = c(-0.1,0.2), breaks = c(-0.1,-0.05,0,0.05,0.1,0.15,0.2), labels = c("-0.1","-0.05","0","0.05","0.1","0.15", "0.2")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  guides(col = guide_legend(nrow = 2)) +
  ggtitle("Switzerland (SHP)")
p_shp <- move_layers(p_shp, "GeomPoint", position = "top")
p_shp <- move_layers(p_shp, "GeomErrorbar", position = "top")
                
p_soep <- dwplot(soep, vline = geom_vline(xintercept=c(0), color = "grey60")) +
  scale_colour_grey(start = .1, end = .1)  + 
  geom_vline(xintercept=0.05, linetype = "dashed", color = "grey60") +
  geom_vline(xintercept=0.1, linetype = "twodash", color = "grey60") +
  scale_x_continuous(limits = c(-0.1,0.2), breaks = c(-0.1,-0.05,0,0.05,0.1,0.15,0.2), labels = c("-0.1","-0.05","0","0.05","0.1","0.15", "0.2")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("Germany (SOEP)")
p_soep <- move_layers(p_soep, "GeomPoint", position = "top")
p_soep <- move_layers(p_soep, "GeomErrorbar", position = "top")

p_us <- dwplot(us, vline = geom_vline(xintercept=c(0), color = "grey60")) +
  scale_colour_grey(start = .1, end = .1)  + 
  geom_vline(xintercept=0.05, linetype = "dashed", color = "grey60") +
  geom_vline(xintercept=0.1, linetype = "twodash", color = "grey60") +
  scale_x_continuous(limits = c(-0.1,0.2), breaks = c(-0.1,-0.05,0,0.05,0.1,0.15,0.2), labels = c("-0.1","-0.05","0","0.05","0.1","0.15", "0.2")) + 
  theme_bw() + 
  theme(legend.position = "none", legend.title=element_blank()) +
  ggtitle("United States")
p_us <- move_layers(p_us, "GeomPoint", position = "top")
p_us <- move_layers(p_us, "GeomErrorbar", position = "top")
p_us


plot_hhinc <- ggarrange(p_shp, p_liss, p_gles, p_soep, p_polat, p_bhps, p_us, nrow = 4, ncol=2, align='v')
ggsave("./Out/Figures/hybrid_eqincdec_ET_col.pdf", plot = plot_hhinc, height = 3 , width = 2.5, scale=4.0)
ggsave("./Out/Figures/hybrid_eqincdec_ET_col.emf", plot = plot_hhinc, height = 3 , width = 2.5, scale=4.0)





